Lab 5 - Data Variation and Co-Variation

Instructions and Overview

In this assignment, we will practice graphically representing variation and covariation in a dataset. To begin you will need to import and clean your dataset. You may reference last week’s lab to help with this. After this, you should follow the prompts and complete the short answer questions. At various points in this assignment, you will be asked to draw insights about your dataset after calling certain functions. When I ask you to draw an insight, I’m not asking you to describe what the function does or to state the results that you get. Instead I’m asking you to interpret those results and consider what this might tell us about the issues represented in the dataset or if it might signal issues of data quality. For instance, stating “the maximum value in the age column is 999,” is not an insight. Instead you should say, “the maximum value in the age column is 999, which is much higher than I would expect and may signal that the data was input wrong or that the data collectors at using 999 to represent null values.”

Getting Started

Load the relevant libraries

library(tidyverse)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
── Attaching packages ─────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.2.0     ✔ purrr   0.3.3
✔ tibble  2.1.3     ✔ dplyr   0.8.3
✔ tidyr   0.8.3     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.4.0
── Conflicts ────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(lubridate)

Attaching package: ‘lubridate’

The following object is masked from ‘package:base’:

    date

library(shiny)
library(shinydashboard)

Attaching package: ‘shinydashboard’

The following object is masked from ‘package:graphics’:

    box

Import and clean example datasets

hospitals <- read.csv("/Users/lpoirier/Documents/GitHub/STS-115/datasets/Hospitals.csv", stringsAsFactors = FALSE)

world_health_econ <- read.csv("https://raw.githubusercontent.com/lindsaypoirier/STS-115/master/datasets/world_health_econ.csv", stringsAsFactors = FALSE)

cases <- read.csv("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv", stringsAsFactors = FALSE)

#Do not worry about this line of code for now. Since the cases data gets appended every day with a new column representing that day's case counts, if we want the total cases per country, we need to add up all of the previous day's counts into a new column. The column below does this for us. 
cases <- 
  cases %>% 
  mutate(Total.Cases = 
           cases %>% 
           select(starts_with("X")) %>% 
           rowSums()
         ) %>%
  select(Province.State, Country.Region, Total.Cases)

hospitals$ZIP <- as.character(hospitals$ZIP)

hospitals$ZIP <- str_pad(hospitals$ZIP, 5, pad = "0") 

is.na(hospitals) <- hospitals == "NOT AVAILABLE"
is.na(hospitals) <- hospitals == -999
is.na(cases) <- cases == ""

hospitals$SOURCEDATE <- ymd_hms(hospitals$SOURCEDATE)
hospitals$VAL_DATE <- ymd_hms(hospitals$VAL_DATE)

Import and clean your dataset.

#Uncomment the lines associated with your datasets below and then run the code chunk.

#dom_violence_calls <- read.csv("https://data-openjustice.doj.ca.gov/sites/default/files/dataset/2019-06/DVRCA_2001-2018.csv", stringsAsFactors = FALSE)
#dom_violence_calls <- dom_violence_calls %>% unite("YEAR_MONTH", YEAR, MONTH, sep = "-", remove = TRUE) %>% mutate(YEAR_MONTH = ymd(YEAR_MONTH, truncated = 1))

#internet_master_plan <- read.csv("https://data.cityofnewyork.us/resource/fg5j-q5nk.csv", stringsAsFactors = FALSE)
#is.na(internet_master_plan) <- internet_master_plan == "N/A"

ipps <- read.csv("https://data.cms.gov/api/views/tcsp-6e99/rows.csv?$limit=200000", stringsAsFactors = FALSE)

#snap <- read.csv("https://github.com/lindsaypoirier/STS-115/blob/master/datasets/SNAP_Policy_Database.csv?raw=true", stringsAsFactors = FALSE)
#is.na(snap) <- snap == -9 
#snap$yearmonth <- ymd(snap$yearmonth, truncated = 1)  

#self_sufficiency_ca <- read.csv("https://github.com/lindsaypoirier/STS-115/blob/master/datasets/CA2020_2018_SSS.csv?raw=true", stringsAsFactors = FALSE)
#self_sufficiency_ca[11:24] <- lapply(self_sufficiency_ca[11:24], gsub, pattern = "[\\$,%() ]", replacement = "")
#self_sufficiency_ca[11:24] <- lapply(self_sufficiency_ca[11:24], as.numeric)

#aqi <- read.csv("https://github.com/lindsaypoirier/STS-115/blob/master/datasets/daily_aqi_by_county_2019.csv?raw=true", stringsAsFactors = FALSE)
#aqi$date <- mdy(aqi$date) 

Exploring Subsetted or Grouped Data

Last week, we explored how your data was defined, how it was categorized, and how missing values were represented. Through this exercise, we were able to develop a general overview of the observations, variables, and values in our datasets. Knowing this we likely have a better sense of areas we want to be able to zoom into and explore in our data. We also sometimes want to be able to zoom out - summarizing something across groups of observations in our data. As we prepare to visualize variation and co-variation in the dataset, let’s first go over how we can both zoom in and zoom out on our datasets.

Zooming in - Exploring Filtered Categorical Data

Filtering to a Category

Filtering is one way we can zoom in on our data - exploring only those observations that meet a particular criteria.

For instance, last week we learned above that one criteria for being designated as a Critical Access Hospital is that the hospital must have 25 or fewer inpatient beds. We may want to see how the values for BEDS change when we filter to (or zoom into) just those observations representing critical access hospitals. Notice how we can combine any number of dplyr verbs using the pipe (%>%). Below I will filter to rows that meet a criteria and then select a variable to examine from the filtered data.

#Run this code chunk.

#df %>% filter(CATEGORICAL_VARIABLE == "VALUE") %>% select(CATEGORICAL_VARIABLE)
hospitals %>% filter(TYPE == "CRITICAL ACCESS") %>% select(NAME, BEDS)

hospitals %>% filter(TYPE == "CRITICAL ACCESS") %>% select(BEDS) %>% summary()

      BEDS      
 Min.   :  3.0  
 1st Qu.: 22.0  
 Median : 25.0  
 Mean   : 27.7  
 3rd Qu.: 25.0  
 Max.   :286.0  
 NA's   :42     

We can see that there are some hospitals in the US that have been designated as critical access hospitals that have more than 25 beds. Since this does not align with the criteria for critical access hospitals that we discovered in our research, it is likely something that we will want to investigate further.

Filtering to Numeric Observations Above or Below a Threshold

We may also want to see which states have hospitals with more than 1200 beds. To do so, we would filter the data to those observations where BEDS is greater than 1200. The first time I do this below, I select() the NAME of the hospital and the STATE. This shows every hospital with more than 1200 beds, along with the STATE in which it is located. However, what if I only wanted to know which states had hospitals with more than 1200 beds? I wouldn’t necessarily want to select() STATE in this case because this would show me the state for every observation in my filtered data. As you can see below, this means that CA would appear 7 times! Instead, I would want to check the distinct() STATES in the filtered data. Check out how I do this below.

#Run this code chunk.

#df %>% filter(NUMBERIC_VARIABLE > VALUE) %>% distinct(CATEGORICAL_VARIABLE)

#Lists name and state of hospitals with more than 1200 beds.
hospitals %>% filter(BEDS > 1200) %>% select(NAME, STATE) 

#Lists state for every hospital with more than 1200 beds in the dataset
hospitals %>% filter(BEDS > 1200) %>% select(STATE) 

#Lists the states wtih more than 1200 beds
hospitals %>% filter(BEDS > 1200) %>% distinct(STATE) 

Zooming out - aka grouping common values and summarizing

At times, we are seeking to get a broader picture of what’s going on in our dataset than provided - grouping observations that share a common value and then performing a calculation to summarize something within each of those groups. In other words, sometimes we want to see our data in aggregate. For instance, I may want to know the total number of hospital beds per state. To calculate this, I would need to group all of the hospital observations by state and then sum the total number of beds in each group. In such cases, we can call group_by() to aggregate the observations with common variable values into groups. Then we will call summarize() to perform a calculation within each of those groups. summarize() takes a set of values and a calculation method and returns a single value. For instance, if we call summarize() with a numeric column in our dataset and “mean” as a calculation method, it will return the average of all the numeric values in that column. When called in conjunction with group_by(), summarize() takes a set of values for each group and a calculation method and returns a single value for each group.

For the hospitals dataset, we will group the observations by STATE and then use summarize to calculate the sum of BEDS per state. Notice below how we are choosing to ignore NA values above by calling na.rm = TRUE. When we do so, we need to keep in mind that we are not summarizing across all observations in the dataset, but only those for which there is a value listed in the variable we are operating on. Because of this, I also calculate the number of observations in each group, the number of observations where the BEDS variable is missing, and the percentage of observations in the group where the variable is missing. This provides important context for how readily we can rely on these numbers. For instance, when you run the code below, note how in Alaska, the number of beds are missing for 28% of the observations.

#Run this code chunk.

#df %>% group_by(CATEGORICAL_VARIABLE) %>% summarize(NEW_VARIABLE_NAME = sum(NUMBERIC_VARIABLE, na.rm = TRUE)) %>% ungroup()

hospitals %>% 
  group_by(STATE) %>% #Group observations by state
  summarize(
    STATES_BEDS = sum(BEDS, na.rm = TRUE), #Calculate the sum of BEDS within each STATE group
    OBSERVATIONS = n(), #Calculate how many observations are in each STATE group
    MISSING_BEDS = sum(is.na(BEDS)), #Calculate how many NAs are in the BEDS variable in each STATE group
    PERCENT_MISSING = sum(is.na(BEDS))/n()*100) %>% #Divide the two values you just calculated to determine the percent of missing data
  ungroup()

Notice that I close each of these calls with ungroup(). When we group_by() a variable, any subsequent function calls will continue to be performed on the grouped data, unless we ungroup() it. This can be important if we want to filter to specific values after we summarize() the data. Assuming that we don’t want to perform a filter operation within each group but on the entire new dataframe created after summarizing, we need to ungroup() the data before performing the filter() operation.

From this function, we see the number of beds across all hospitals per state. Depending on the question we are asking, this may or may not be relevant. For instance, if I’m wondering how much hospital infrastructure is available to support Covid-19 patients, one (of a number of factors) I need to consider before presenting this data is which types of hospitals are accepting Covid-19 patients. Are rehabilitation hospitals accepting patients? Psychiatric hospitals? Military hospitals? If they aren’t now, will they at some point? Further, some states are talking about cordoning off hotels for Covid-19 patients. How do we account for this change in the number of hospital beds (something definitely not represented in our data based on the way hospital has been defined). We need to do external research to answer these questions. Then we may wish to filter our data to relevant hospital types. For instance, at this moment, we may filter our data to only include beds at General Acute Care Hospitals. We also know that some hospitals in the dataset are closed. We need to also filter these out before presenting the data. Notice how below, I can do this by simply copying and pasting the code from above and adding one filter statement before grouping the data.

#Run this code chunk.

hospitals %>% 
  filter(TYPE == "GENERAL ACUTE CARE" & STATUS == "OPEN") %>%
  group_by(STATE) %>% #Group observations by state
  summarize(
    STATES_BEDS = sum(BEDS, na.rm = TRUE), #Calculate the sum of BEDS within each STATE group
    OBSERVATIONS = n(), #Calculate how many observations are in each STATE group
    MISSING_BEDS = sum(is.na(BEDS)), #Calculate how many NAs are in the BEDS variable in each STATE group
    PERCENT_MISSING = sum(is.na(BEDS))/n()*100) %>% #Divide the two values you just calculated to determine the percent of missing data
  ungroup()

In other words, often times to answer questions within a dataset, we need to both zoom in and out on data - honing in on certain observations and then generalizing across them. We cannot answer questions well if we don’t have a good understanding of what’s included in our data and how issues are defined. Had we not known that hospitals that are closed and hospitals that are classed as rehabs or psychiatric facilities were included in the data, we may have made some poor assumptions about the number of beds available. Also note how, in every step of data analysis, we have to make decisions about what to include and what to exclude in the analysis. Data analysts play a very active role in shaping the knowledge that gets produced from data. The numbers can never speak for themselves.

When I Need to Zoom In or Out

For some of you, these functions will be necessary to employ before performing operations across numeric variables in your dataset. This is because, as we learned last week, some of you have observational units that span multiple time periods, multiple geographies, or multiple issues. Before performing an operation across a numeric variable, we need to ensure all of the values in that variable are referring to observations reported across the same timeframe or geographic scale.

In the hospitals dataset, this is less of a concern because as we learned last week, the observational unit of the dataset is one thing without qualifiers - a hospital. The BEDS variable is always going to refer to the number of BEDS at a hospital.

With the world_health_econ data, we learned last week that every observation refers to a country and a year. Let’s say that we wanted to call summary() on life_exp variable to compare life expectancies across countries. Without first zooming in to to a specific year, we would be including multiple values taken at the same place at different times. Let’s filter the data to only include the most recent reporting year and then call summary():

#Run this code chunk.

world_health_econ %>%
  filter(year == max(year, na.rm = TRUE)) %>% #Note that this is how we can fliter to the rows with the maximum value in a variable; in this case, this would be the most recent year. 
  select(life_exp) %>%
  summary()

    life_exp    
 Min.   :32.50  
 1st Qu.:63.95  
 Median :73.35  
 Mean   :70.52  
 3rd Qu.:77.05  
 Max.   :83.20  
 NA's   :6      

In rarer cases, we will need to group_by() and summarize() our data before statistically analyzing it. This is the case with the cases dataset. Each observation in the cases dataset refers to a combination of things - a province and a country. The issue is that not every country is reporting data at the province level. This means that sometimes the Total.Cases variable is referring to the number of cases in a country and sometimes it is referring to the number of cases in a Province. If we wanted to statistically analyze the number of cases across observations, we first need to transform the dataset so that each observation is reported at the same scale. In other words, we need to zoom out to standardize the observational unit across the dataset to the country level (since we do not always have data at the province level). To do that, we will need to group_by() Country.Region and then summarize the sum of Total.Cases.

#Run this code chunk.

cases %>%
  group_by(Country.Region) %>%
  summarize(Total.Cases = sum(Total.Cases, na.rm = TRUE)) %>%
  ungroup()

Once we’ve done this, we can confidently statistically analyze the Total.Cases variable because we know that the number is always referring to the number of cases in a country.

Alternatively, we could also filter to one country to analyze numeric variables across provinces in that country. If I were to call summary() on the Total.Cases variable in the cases dataset, I would be gathering statistics across numbers reported at different geographic scales. Below, I filter the cases dataset to one Country - Canada - so that I’m gathering statistics across all provinces in Canada.

#Run this code chunk.

cases %>% 
  filter(Country.Region == "Canada") %>%
  select(Total.Cases) %>%
  summary()

  Total.Cases     
 Min.   :     16  
 1st Qu.:    906  
 Median :  12245  
 Mean   : 147491  
 3rd Qu.:  78285  
 Max.   :1073165  

Let’s test the extent to which you will need to zoom in or zoom out to statistically analyze numeric values your data. I’ve selected a numeric variable for each of your datasets. Comment out the line with your dataset and complete the statement below. Think what variables make up the unique key in your datasets. If you have more than one variable in your unique key, make sure that each is represented in your statement below. As a hint, I will note that only one of the six groups has only one variable in their unique key ((ahem infrastructure). Three of the groups have two variables, and two groups (ahem income security) have three variables.

Environment group only: On import your data does require three variables to constitute a unique key. However, with the code below, I’m going to create a new variable in your dataset to uniquely refer to each county so that you do not have reference both the state code and the county code to refer to a county. This variable will be called fips, and it will appear as the last variable in the dataset.

#Environment group only: uncomment this line to create this new variable. Respond to the questions below with this change in the dataset.

#aqi <- aqi %>% mutate(fips = paste(str_pad(as.character(stateCode), 2, side = "left", "0"), str_pad(as.character(countyCode), 3, side = "left", "0"), sep =""))

#Uncomment the line associated with your dataset and fill in the blank. Then run the code.

#paste(df$NUMERIC_VARIABLE[1], "refers to a number/measure of [FILL NUMERIC VARIABLE] in a _____ in my dataset.")

#Example:
paste(hospitals$BEDS[1], "refers to a number of beds in a hospital in my dataset.")

[1] "25 refers to a number of beds in a hospital in my dataset."

paste(world_health_econ$pop[1], "refers to a number of people in a country in a given year in my dataset.")

[1] "18100000 refers to a number of people in a country in a given year in my dataset."

paste(cases$Total.Cases[1], "refers to a number of cases in a country and/or province in my dataset.")

[1] "94226 refers to a number of cases in a country and/or province in my dataset."

#paste(dom_violence_calls$TOTAL_CALLS[1], "refers to a number of calls in a _____ in my dataset.")
#paste(aqi$numberSitesReporting[1], "refers to the number of sites reporting in a _____ in my dataset.")
#paste(internet_master_plan$total_number_of_households[1], "refers to a number of households in a _____ in my dataset.")
paste(ipps$Total.Discharges[1], "refers to a number of discharges in a given hospital based on a diagnosis classification in my dataset.")

[1] "28 refers to a number of discharges in a given hospital based on a diagnosis classification in my dataset."

#paste(snap$certearnavg[1], "refers to a measure of earnings in a _____ in my dataset.")
#paste(self_sufficiency_ca$housingCosts[1], "refers to a measure of housing costs in a _____ in my dataset.")

How did you fill in the last _____?

Is your observational unit one thing (e.g. one hospital, or one country)? If this is the case, it will likely not be as essential for you to zoom in or zoom out before operating on numeric variables. If so, we will say that you do not have qualified units of observation.

OR

Is your observational unit a combination of things or factors (e.g. one chemical reported at a particular facility or one census tract reporting in a particular year)? If this is the case, it will likely be essential for you to zoom in or zoom out before operating on numeric variables. If so, we will say that you do have qualified units of observation.

For each of the groups that has qualified units of observation (most of you do), in some places throughout this lab, you will need to filter your data to particular observations before analyzing across a variable. This is because we will be summarizing information across groups of data, and it will be important to ensure that you are summarizing information across like observations.

Think about what you might filter to in order to ensure that you will be comparing like observations (hint: it will involve a variable in your unique key). Perhaps you will filter to the most recent year, so that you can compare observations across geographies in that year. Or perhaps you will filter to a particular geography, so that you can compare observations across time in that geography. Or perhaps you will filter to a particular diagnosis group, so that you can compare costs across hospitals for that diagnosis. Or perhaps you will filter to a particular year and family type so that you can compare observations across counties in that year for that family type. Characterize one way you might filter your data below. Be specific. Which variable in the dataset will you filter on and to what value(s) will you filter it to?

I will filter my data to the variables DRG.Definition, and I will filter the values to chest pain diagnoses across hospitals.

Zooming in and/or Zooming Out on Your Own Data

Select one of the values that you identified from calling distinct() on a categorical variable in last week’s lab. Filter the dataset to the rows representing that value, select a numeric variable to explore, and then call summary(). If you have qualified units of observation, be sure to first zoom into a set of observations in your data (using filter()).

#Uncomment the appropriate lines below, and fill in your data frame, variables, and value. Run the code.

#_____ %>% filter(_____ == "_____") %>% select(_____) %>% summary()

#ipps %>% filter(DRG.Definition == "313 - CHEST PAIN") %>% select(Total.Discharges) %>% summary()

#If you have qualified units of observation
ipps %>% filter(DRG.Definition == "313 - CHEST PAIN") %>% select(Average.Total.Payments) %>% summary()

 Average.Total.Payments
 Min.   : 3392         
 1st Qu.: 4326         
 Median : 4940         
 Mean   : 5431         
 3rd Qu.: 5994         
 Max.   :14467         

#Income security group only
#_____ %>% filter(_____ == _____ & _____ == _____ & _____ == "_____") %>% select(_____) %>% summary()

What question might this analysis help to address?

A question that this analysis might help me to address would be what is the variation in the average total payments for chest pain related diagnosis across hospitals in the US.

Are there any other variables in your dataset that you need to take into consideration before directing this analysis towards answering that question? In other words, do you need to zoom into any specific areas of the dataset (by filtering) in order to appropriately address this question? If so, which? (For instance, for the hospitals dataset we needed to filter to hospitals that were open in order to address questions about infrastructure). Be sure to adjust your code above to reflect this.

No other variables in my dataset need to be considered.

What else would we need to know to fully address this question? Here you may consider what you know about how this dataset was produced and its limitations.

There are limitations to this dataset. For instance, we would need to know exactly what conditions would be classified as "Chest Pain", and it would be helpful to know why this was the primary diagnosis and whether there were other symptoms associated with it as well. Furthermore, it would be helpful to know what tratments are done for this type of diagnosis, and whether each hospital has different criteria when diagnosing chest pain.

What insight can you draw from calling summary() on your filtered dataset?

One insight that I can draw from calling summary() is that there is a significant difference between the minimum and maximum average payments for chest pain disgnosis. This might be due to different treatments performed, lengths of stay of patients, and variation in hospital expenses.

Select a numeric variable in your dataset that represents the extent or scale of the issue you are studying. Pick a number that you believe serves as a good indicator that this issue is at a notable extent or scale, and filter the dataset to all the rows greater than (or less than) this number. Check the remaining distinct values in a categorical variable in the dataset. If you have qualified units of observation, be sure to first zoom into a set of observations in your data (using filter()).

#Uncomment the appropriate lines below, and fill in your data frame, variables, condition, and value. Run the code.
#_____ %>% filter(_____ _____ _____) %>% distinct(_____)

#If you have qualified units of observation
ipps %>% filter(DRG.Definition == "313 - CHEST PAIN" & Total.Discharges > 180) %>% distinct(Provider.Name, Total.Discharges, Provider.City, Provider.State)

#Income security group only
#_____ %>% filter(_____ == _____ & _____ == _____ & _____ _____ _____) %>% distinct(_____)

What question might this analysis help to address?

This analysis might help me address what are the hospitals with the highest number of paients dismissed that had a chest pain diagnosis. 

Are there any other variables in your dataset that you need to take into consideration before directing this analysis towards answering that question? In other words, do you need to zoom into any specific areas of the dataset (by filtering) in order to appropriately address this question? If so, which? Be sure to adjust your code above to reflect this.

It would be helpful to know where these hospitals are located and in what US State.  

What else would we need to know to fully address this question? Here you may consider what you know about how this dataset was produced and its limitations.

It would be helpful to see how these hospitals classify chest pain diagnosis, and what type of criteria they use.

What insight can you draw from calling distinct on the filtered data?

From calling distinct on my filtered data I can see the name and city of the few hospitals that have more than 180 discharges of patients diagnosed with chest pain, and I can also see that they are located in New York, Texas, and Michigan. This could help me narrow down my research to try explain why these hopsitals have such a high number of chest pain diagnoses, and whether it has to do with the way they classify it.

Select a categorical variable that you would like to group your data by, so that you can summarize some statistics across each grouping. You may group your data by a particular year, by a particular location (such as a state or a region), or by a particular category. Then select a numeric variable in your dataset to summarize by. For instance, you may want to sum the total number of reports in a given year, or find the average number of cases reported in a certain state. Also calculate the percent of observations where data is missing in each group. If you have qualified units of observation, you may want to group the data by one of the variables in your unique key. For instance, if your unique key is a county and year, then perhaps you want to group the data by county and summarize something across each year. If you have qualified units of observation and choose to group by a variable that is not in your unique key, then be sure to filter the data as you have been above.

#Uncomment the appropriate lines below, and fill in your data frame, variables, and summarize variable name, and math function. Then run the code.

ipps %>% group_by(DRG.Definition) %>% summarize(Total.Discharges.DRG = sum(Total.Discharges, na.rm = TRUE), OBSERVATIONS = n(), MISSING = sum(is.na(Total.Discharges)), PERCENT_MISSING = sum(is.na(Total.Discharges))/n()*100)

#If you have qualified units of observation (and not grouping by the qualifier).
#_____ %>% filter(_____ == _____) %>% group_by(_____) %>% summarize(_____ = _____(_____, na.rm = TRUE), OBSERVATIONS = n(), MISSING = sum(is.na(_____)), PERCENT_MISSING = sum(is.na(BEDS))/n()*100)

What question might this analysis help to address?

This analysis could mainly helping me address the question about how many types of classifications diagnoses exist in this dataset, and what the total number of patients diagnosed with each condition is. 

Are there any other variables in your dataset that you need to take into consideration before directing this analysis towards answering that question? In other words, do you need to zoom into any specific areas of the dataset (by filtering) in order to appropriately address this question? If so, which? Be sure to adjust your plot above to reflect this.

No, there are not.

What else would we need to know to fully address this question? Here you may consider what you know about how this dataset was produced and its limitations.

I would need to acknowledge that there are so many possible variations in classifying diagnoses, and that many of these classifications potentially group together different conditions in one diagnosis. If this dataset were to consider all possible variations and conditions, the number of diasgnoses would much much higher. Also, patients that have been classified with one of these diagnoses might also have other conditions that are not registered, as patients only have a primary diagnosis registered. 

What insight can you draw from grouping and summarizing?

I can see what the most common diagnoses are, and how many patients receive that diagnosis. Also, I can see what the rarest diagnoses and procedures are. Also, I can see from the number of observations how many hospitals have at least 1 specific diagnosis. For instance, I can see that there are a total of 28 patients diagnosed with AUTOLOGOUS BONE MARROW TRANSPLANT W/O CC/MCC, but only 2 observations, which means that only 2 hospitals have performed this type of procedure. This could tell me that this type of procedure associated with this diagnosis is pretty rare and might be only performed in a couple of hospitals in the US which might be specialized and equipped to treat such condition. 

Combine any combination of the 4 verbs we learned in class this week or last (select, filter, group by, or summarize) to explore your dataset further. You may also use arrange, summary, or distinct.

ipps %>% filter(DRG.Definition == "017 - AUTOLOGOUS BONE MARROW TRANSPLANT W/O CC/MCC") %>% distinct(Provider.Name, Total.Discharges, Provider.City, Provider.State)

NA

What question might this analysis help to address?

What hospitals are specialized have treated AUTOLOGOUS BONE MARROW TRANSPLANT W/O CC/MCC.

Are there any other variables in your dataset that you need to take into consideration before directing this analysis towards answering that question? In other words, do you need to zoom into any specific areas of the dataset (by filtering) in order to appropriately address this question? If so, which? Be sure to adjust your plot above to reflect this.

No. 

What else would we need to know to fully address this question? Here you may consider what you know about how this dataset was produced and its limitations.

Again, we would need to keep into consideration the fact that these patients might have had other conditions as well, but this was the primary diagnosis and treatment. 

What insight can you draw from running this function?

I could potentially deduct that these two hospitals are specialized in treating this condition, although more research should be done to find out more about the type of procedures and treatments available at other hospitals as well.  

---

Variation

Variation is the extent to which the values in a particular variable vary from observation to observation. Examining variation involves looking at the distribution of values in a particular column in the dataset. Do we have a whole bunch of one particular value in a certain variable, and very few of another? Or maybe, do we have a more even distribution of values across a variable?

ggplot

At this point in the assignment, we will begin leveraging the Tidyverse package ggplot to create plots for visualizing the data. To create a plot with ggplot, we will follow this basic formula:

df %>% 
  ggplot(aes(x = VARIABLE_NAME)) + 
  CHART_TYPE

For example, for a bar chart, you will call:

df %>% 
  ggplot(aes(x = VARIABLE_NAME)) + 
  geom_bar()

For a column chart, you will call:

df %>% 
  ggplot(aes(x = VARIABLE_NAME, y = VARIABLE_NAME)) + 
  geom_col()

Let’s break that down a bit. First, you will call your dataframe. Following your dataframe and a pipe, you will call ggplot(), which basically tells R to prepare to create a plot. Inside ggplot, you will list aesthetics. These are variables in the dataset that you would like to appear on your plot. Setting x = VARIABLE_NAME tells ggplot() what variable to plot on the x-axis. Setting y = VARIABLE_NAME tells ggplot() what variable to plot on the y-axis. Finally, following a plus (+) sign, you tell ggplot which type of plot to create. The ggplot cheatsheet lists a number of plots that you can create with ggplot, as well as a number of different ways to style the plot. We will practice several of these below.

For every plot that you produce, I will expect you to add a title and labels to the x and y axis. You can do this as follows:

df %>% 
  ggplot(aes(x = VARIABLE_NAME, y = VARIABLE_NAME)) + 
  geom_col() +
  labs(title = "FILL TITLE", x = "FILL X-AXIS LABEL", y = "FILL Y-AXIS LABEL)

There are also a number of useful tools for styling your plots. For instance we can set the theme of the plot to look a bit more polished by adding “+ theme_bw()” to the plot. I will do this for all plots in this lab.

df %>% 
  ggplot(aes(x = VARIABLE_NAME, y = VARIABLE_NAME)) + 
  geom_col() +
  labs(title = "FILL TITLE", x = "FILL X-AXIS LABEL", y = "FILL Y-AXIS LABEL") +
  theme_bw()

Two styling issues that I’m sure will come up in most of your plots include: * changing x or y axis tick numbers from scientific to comma notation: + scale_x_continuous(labels = scales::comma) OR + scale_y_continuous(labels = scales::comma) * turning x axis tick marks 90 degrees so that they do not overlap: + theme(axis.text.x = element_text(angle = 90, hjust=1)) OR + theme(axis.text.y = element_text(angle = 90, hjust=1))

df %>% 
  ggplot(aes(x = VARIABLE_NAME, y = VARIABLE_NAME)) + 
  geom_col() +
  labs(title = "FILL TITLE", x = "FILL X-AXIS LABEL", y = "FILL Y-AXIS LABEL") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust=1)) + #Turn labels 90 degrees
  scale_x_continuous(labels = scales::comma) #Change labels from scientific to comma notation

So how do we visualize variation with ggplot? Below I describe two different plots that you can leverage to visualize variation - a bar plot and a frequency plot.

Bar Plot

A bar plot displays the number of times each value appears in a categorical variable. This will tell us how the observations in the dataset vary in regards to that variable. In other words, this plot will communicate the number of observations in your dataset by that variable.

#Run this code chunk.

#df %>% ggplot(aes(x = CATEGORICAL_VARIABLE)) + geom_bar() + labs(title = "TITLE", x = "X-AXIS NAME", y = "Y-AXIS NAME")

hospitals %>% 
  ggplot(aes(x = TYPE)) + 
  geom_bar() +
  labs(title = "Number of Hospitals in the US by Type", x = "Type", y = "Count of Hospitals") + #Adds a title to the plot
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust=1)) #Changes x-axis tick labels 90 degrees 

Remember that this dataset includes hospitals that are designated as closed. Depending on the question we are trying to address, we may wish to zoom in to only the observations signifying a hospital that is open before creating this plot.

#Run this code chunk.

#df %>% ggplot(aes(x = CATEGORICAL_VARIABLE)) + geom_bar() + labs(title = "TITLE", x = "X-AXIS NAME", y = "Y-AXIS NAME")

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = TYPE)) + 
  geom_bar() +
  labs(title = "Number of Hospitals in the US that are Open by Type", x = "Type", y = "Count of Hospitals") + #Adds a title to the plot
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust=1)) #Changes x-axis tick labels 90 degrees 

Titling a Bar Plot

Note how I titled my first plot above: “Number of Hospitals in the US by Type” Remember last week, when we identified what makes each observation in our dataset unique? Here I am counting the observations by Type, and in order to know what I’m counting, I need to know what each observation refers to. A good formula for titling bar plots is as follows:

Number of _____ by [x-axis variable name]

In order to fill in the blank line above, consider the statement we produced last week: “I have nrow(df) unique _____ represented in my dataset.” That blank line told us what each observation in the dataset represented - or its observational unit. However we filled in that blank line should also be how we fill in the title of a bar plot.

The [x-axis variable name] should be your x-label and “Count of ______” (filled the same as above) should be your y-label. Note that if you filter your dataset, you should account for this in the title: “Number of Hospitals in the US by Type”

What if I have qualified units of observation?

When we create a bar plot, we are counting the number of observations that fall into each category. If there is only one variable that makes up your unique key, that one variable will represent what is being counted (e.g. above we are counting hospitals by category). However, if there are multiple variables in your unique key, then identifying what it is that you are counting becomes a little more complicated. For instance, let’s say your data reports the population of each country each year as it does in the world_health_econ dataset. Now let’s say that you wanted to plot the number of countries per continent. If you were to call:

#Run this code chunk.

world_health_econ %>% 
  ggplot(aes(x = continent)) + 
  geom_bar() +
  labs(title = "Number of Countries per Continent - Incorrect", x = "Continent", y = "Count of Countries") + #Adds a title to the plot
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust=1)) #Changes x-axis tick labels 90 degrees 

… you would be counting the combination of the number of countries and years per continent. There are not 800 countries in Africa, but Africa appears in the dataset 800 times because each country in Africa has several rows in the dataset - one for each reporting year. In other words, each country is represented in the bar for every year that it was included in the dataset. The y-axis does not just represent countries but both countries and years. If we want the y-axis to only be counting one thing, then we need to first reduce the dataset to values in a particular context. You can do this by filtering the data as you had been doing above.

#Run this code chunk.

world_health_econ %>% 
  filter(year == max(year, na.rm = TRUE)) %>%
  ggplot(aes(x = continent)) + 
  geom_bar() +
  labs(title = "Number of Countries per Continent in the most Recent Reporting Year", x = "Continent", y = "Count of Countries") + #Adds a title to the plot
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust=1)) #Changes x-axis tick labels 90 degrees

Note the addition to my title above. If you filter your dataset, the formula for titling changes a bit to Frequency of _____ across [x-axis variable name] in [filtered value]

Select a categorical variable for which you want to visualize the frequency of times it appears in the dataset.

I recommend that you select one of the same categorical variables that you analyzed with the distinct() function last week. If you have qualified units of observation, be sure to filter your data first so that all of the observations you are counting have one variables as a unique key. If you notice styling issues with your data, be sure to consider my notes about how to fix them above.

#Uncomment the line below and fill appropriately. Add a title and labels to your plots, and adjust its style to be legible. Then run the code.
#_____ %>% ggplot(aes(x = _____)) + geom_bar()

#If qualified:
ipps %>% filter(DRG.Definition == "313 - CHEST PAIN") %>% ggplot(aes(x = Provider.State)) + geom_bar() +
labs(title = "Number of Chest Pain Diagnoses across Hospitals by State", x = "State", y = "Count of Chest Pain Diagnoses") +
theme_bw()

Reflect on the distribution of categories in the dataset. Is there an even distribution of observations across each category, or are certain categories more represented than others? Why might this be? What does this say about the social, political, or economic landscape of your topic?

No, there is not an even distribution of observations across States. As it can be seen from the graph there are some States that have a significantly larger number of patients diagnosed and treated for chest pain. One of the reasons why is that some States have higher number of people that are in Medicare.

Reflect on what you learned about the history of the categories you plotted above in last week’s lab. How have the social, political, and economic forces shaping this categorization impacted how we count observations related to this topic? How might this plot look different if this variable had been categorized in a different way?

Fill response here. 

Frequency Plot

A frequency plot will display the distribution of values in a numeric variable within a designated set of increments. This will tell us how the observations in the dataset vary in regards to that variable. In other words, this plot will communicate the number of observations in your dataset broken down into increments of the values in that variable.

Consider the plot below. This plot tells us the distribution of beds across open hospitals. It tells us how many hospitals there are in each increment of 10 beds.

#Run this code chunk.

#df %>% ggplot(aes(x = NUMERIC_VARIABLE)) + geom_freqpoly(binwidth = 1) + labs(title = "TITLE", x = "X-AXIS NAME", y = "Y-AXIS NAME")

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = BEDS)) +
  geom_freqpoly(binwidth = 10) +
  labs(title = "Distribution of Beds across Hospitals in the US that are Open", x = "Beds", y = "Count of Hospitals") +
  theme_bw()

Note that binwidth refers to the size of the increments at which frequency will be calculated. Above the binwidth is set to 10. This means that ggplot will display the frequency of each value at intervals of 10, 20, 30, 40, etc. When we set the bindwidth to 1, ggplot will display the frequency of each value at intervals of 1, 2, 3, 4. etc. What difference does this make?

Notice what happens when we set the binwidth to 1. While above we count the number of hospitals with 0-10 beds 10-20 beds, 20-30 beds, etc, this will count the number of hospitals with 0-1 beds, 1-2 beds, 2-3 beds, and so on. Because we are counting the number in such small increments, the plot will look much more jagged and will take a longer time to load. This plot displays the counts in finer granularity than the first plot.

#Run this code chunk.

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = BEDS)) +
  geom_freqpoly(binwidth = 1) +
  labs(title = "Distribution of Beds across Hospitals in the US that are Open", x = "Beds", y = "Count of Hospitals") +
  theme_bw()

When we set the binwidth to 100, we count the number of hospitals with 0-100 beds, 100-200 beds, 200-300 beds, and so on. Because we are counting the number in larger increments, the plot will look much smoother and will take less time to load. This plot displays the counts in thicker granularity than the first plot.

#Run this code chunk.

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = BEDS)) +
  geom_freqpoly(binwidth = 100) +
  labs(title = "Distribution of Beds across Hospitals in the US that are Open", x = "Beds", y = "Count of Hospitals") +
  theme_bw()

Titling a Frequency Plot

Note how I titled my plot above: “Distribution of Beds across Hospitals in the US that are Open” Consider again what makes each observation unique. A good formula for titling frequency plots is as follows:

Frequency of [x-axis variable name] across _____

Again, we can fill in the blank with our observational unit. The [x-axis variable name] should be your x-label and “Count of ______” (filled the same as above) should be your y-label.

What if I have qualified units of observation?

Let’s talk about what would happen if I were to make a frequency plot of the Total.Cases in the cases dataset:

#Run this code chunk.

cases %>% 
  ggplot(aes(x = Total.Cases)) + 
  geom_freqpoly(binwidth = 10000) +
  labs(title = "Distribution of Cases across ____", x = "Total Cases", y = "Count of _____") + # To add titles and labels
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust=1)) + #Turn labels 90 degrees
  scale_x_continuous(labels = scales::comma) #Change labels from scientific to comma notation

We know that each observation in the cases dataset refers to a province/country pair. So here I’m counting the number of province/countries in each bracket of total cases. Since only some countries in this dataset are broken into provinces, we are comparing counts of cases across different geographic scales - sometimes at the province level and sometimes at the country level. In this case, it makes more sense to total up the number of cases per country and then plot the distribution of cases across countries. We can use group_by() and summarize() to do this:

#Run this code chunk.

cases %>% 
  group_by(Country.Region) %>%
  summarize(Total.Cases = sum(Total.Cases, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = Total.Cases)) + 
  geom_freqpoly(binwidth = 100000) +
  labs(title = "Distribution of Cases across Countries", x = "Total Cases", y = "Count of Countries") + # To add titles and labels
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust=1)) + #Turn labels 90 degrees
  scale_x_continuous(labels = scales::comma) #Change labels from scientific to comma notation

What about the world_health_econ data? Without separating out these units of observation, we would be visualizing multiple values reported at the same place at different periods in time (i.e. every year since 1995). Instead, we want to zoom into a single year in the dataset so we are just comparing values across place.

#Run this code chunk.

world_health_econ %>% 
  filter(year == max(year, na.rm = TRUE)) %>%
  ggplot(aes(x = tot_health_sp_pp)) + 
  geom_freqpoly(binwidth = 100) +
  labs(title = "Distribution of Total Health Spending per Person across Countries in 2010", x = "Total Health Spending per Person", y = "Count of Countries") + # To add titles and labels
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust=1)) + #Turn labels 90 degrees
  scale_x_continuous(labels = scales::comma) #Change labels from scientific to comma notation

Note the addition to my title above. If you filter your dataset, the formula for titling changes a bit to Frequency of [x-axis variable name] across _____ in [filtered value]

Select a numeric variable for which you want to visualize the distribution of a set of values.

Be sure to select a variable that describes something about the observational unit and not another categorical variable in your dataset. For instance, let’s say you had the following data table - with each row reporting an environmental violation at a facility:

Violation Number | Facility Name | Facility Town | Population of Facility Town _________________ | _________________ | _________________ | _________________ 1234567 | Facility A | Tarrytown | 90000 2345678 | Facility B | Tarrytown | 90000 3456789 | Facility C | Another Town | 70000

In this table, we would not want to create a frequency plot with population of facility town. This is because our observational unit is a violation, not a town, and population does not describe something about the violation but instead describes something about the town the facility is in. If we were to create a frequency plot with this variable, we would be counting the number of times each population value appears in the dataset - something that does not make much sense, since we can have the same town’s population town appear several times in the dataset if there are more than one violations or more than one facilities in a town. None of you should have to worry about this in your datasets. All of your measures describe the observational unit.

If you have qualified units of observation, be sure to first zoom into a set of observations in your data (using filter()).

#Uncomment the line below and fill appropriately. Add a title and labels to your plots. Run the code.
ipps %>% filter(DRG.Definition == "313 - CHEST PAIN") %>% ggplot(aes(x = Total.Discharges))  + geom_freqpoly(binwidth = 10) 

This gives us information about the distribution of values in the dataset. Reflect on the distribution of values. What are the range of values represented in the data? Are the values evenly distributed, or are certain values more represented than others? Why might this be? Do any of the values surprise you? Why?

Fill response here. 

Why did you select the binwidth that you did? How might the story your plot tells change if you were to change the binwidth? What anomalies might be hidden with a larger binwidth, and what trends might be hidden with a smaller binwidth?

Fill response here. 

Check how the numeric variable was defined in the data dictionary, and quote the definition below. How might the counts represented in your frequency plot appear differently if this variable was defined differently?

Fill response here. 

Co-variation

Co-variation is the extent to which the values that constitute two or more variables vary in relation to one another. To visualize co-variation, we might create:

Count Plots

Count plots display how many times two categorical values appear together in a dataset. For instance, in the count plot below, we display the number of open hospitals with each combination of OWNER and TYPE.

#Run this code chunk.

#df %>% ggplot(aes(x = CATEGORICAL_VARIABLE, y = CATEGORICAL_VARIABLE)) + geom_count()

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = TYPE, y = OWNER)) + 
  geom_count() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  labs(title = "Number of Hospitals in the US that are Open by Type and Ownership", x = "Type", y = "Ownership")

Titling a Count Plot

Note how I titled my plot above: “Number of Hospitals in the US that are Open by Type and Ownership” Consider again what makes each observation unique. Here I am counting the observations by Type and Ownership, and in order to know what I’m counting, I need to know what each observation refers to. A good formula for titling count plots is as follows:

Number of _____ by [x-axis variable name] and [y-axis variable name]

The blank should be filled with your unit of observation. The [x-axis variable name] should be your x-label and [y-axis variable name] should be your y-label.

What if I have qualified units of observation?

If this is the case I would encourage you to include one of the qualified variables in the x or y-axis. For instance, if the unique key for world_health_econ is country and year, I can include year as the y-axis below to visualize how counts of observations change over time. (Notice how they don’t in the plot below.) Alternatively, if you wish to compare two categorical variables that are not a part of the unique key, be sure to filter the data so that only one variable constitutes the unique key.

#Run this code chunk.

world_health_econ %>% 
  ggplot(aes(x = continent, y = as.factor(year))) + 
  geom_count() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  labs(title = "Number of Countries by Continent and Year", x = "Continent", y = "Year")

Stacked Frequency Plots

Stacked frequency plots display the distribution of numeric values in a variable, grouped by a categorical value. For instance, the plot below display the distribution of beds across open hospitals, categorized by the hospital type.

#Run this code chunk.

#df %>% ggplot(aes(x = NUMERIC_VARIABLE, col = CATEGORICAL_VARIABLE)) + geom_freqpoly(binwidth = 1)

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = BEDS, col = TYPE)) + 
  geom_freqpoly(binwidth = 100) +
  labs(title = "Frequency of Beds across Hospitals in the US that are Open by Hospital Type", x = "Beds", y = "Count of Hospital", col = "Type") + # To add titles and labels
  theme_bw() 

Stacked frequency plots work best when categorizing by a variable with 10 or fewer distinct values. Otherwise, it can be tricky to see the differences in color gradations. If all of your categorical variables have more than 10 distinct values, one thing you might consider is first filtering your data to a few representative categories. For instance, let’s say that I would like to see the distribution of beds in hospitals across states. Since there are 57 states, if I were to categorize the distribution by state, the plot would be very difficult to read.

#Run this code chunk.

#df %>% ggplot(aes(x = NUMERIC_VARIABLE, col = CATEGORICAL_VARIABLE)) + geom_freqpoly(binwidth = 1)

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = BEDS, col = STATE)) + 
  geom_freqpoly(binwidth = 100) +
  labs(title = "Frequency of Beds across Hospitals in the US that are Open by Hospital State", x = "Beds", y = "Count of Hospital", col = "Type") + # To add titles and labels
  theme_bw() 

In this case, we may wish to first filter our data to a few representative states using %in%.

#Run this code chunk.

#df %>% ggplot(aes(x = NUMERIC_VARIABLE, col = CATEGORICAL_VARIABLE)) + geom_freqpoly(binwidth = 1)

hospitals %>% 
  filter(STATUS == "OPEN" & STATE %in% (c("CA", "MA", "NY", "FL", "LA"))) %>%
  ggplot(aes(x = BEDS, col = STATE)) + 
  geom_freqpoly(binwidth = 100) +
  labs(title = "Frequency of Beds across Hospitals in the US that are Open by Hospital State", x = "Beds", y = "Count of Hospital", col = "Type") + # To add titles and labels
  theme_bw() 

Titling a Stacked Frequency Plot

Note how I titled my plot above: “Frequency of Beds across Hospitals in the US that are Open by Hospital Type” Consider again what makes each observation unique. Here I am counting the observations by number of Beds and Hospital Type. A good formula for titling stacked frequency plots is as follows:

Frequency of [x-axis variable name] across _____ by [col variable name]

The blank should be filled with your unit of observation. The [x-axis variable name] should be your x-label, “Count of ______” (filled the same as above) should be your y-label, and the [col variable name] should be your col-label.

What if I have qualified units of observation?

If this is the case I would encourage you to include one of the qualified variables in the col variable. For instance, if world_health_econ is qualified by country and year, I can include year as the col variable below to visualize how the frequency of countries with various life expectancies changes over time.

#Run this code chunk.

#df %>% ggplot(aes(x = NUMERIC_VARIABLE, col = CATEGORICAL_VARIABLE)) + geom_freqpoly(binwidth = 1)

world_health_econ %>% 
  ggplot(aes(x = life_exp, col = as.factor(year))) + 
  geom_freqpoly(binwidth = 5) +
  theme_bw() +
  labs(title = "Frequency of Life Expectancies across Countries by Year", x = "Life Expectancy", y = "Count of Countries", col = "Life Expectancy") 

Note that this is one plot that is particularly susceptible to missing data. If certain countries did not report data in certain years, the count of countries in a bracket will appear lower, not necessarily because fewer countries fell within a certain life expectancy bracket, but because fewer countries reported that life expectancy.

Alternatively, you could zoom in to one value in one of your qualified variables, filtering to a specific subset of observations and then divide by a different categorical variable.

#Run this code chunk.

#df %>% ggplot(aes(x = NUMERIC_VARIABLE, col = CATEGORICAL_VARIABLE)) + geom_freqpoly(binwidth = 1)

world_health_econ %>% 
  filter(year == max(year, na.rm = TRUE)) %>%
  ggplot(aes(x = life_exp, col = continent)) + 
  geom_freqpoly(binwidth = 5) +
  theme_bw() +
  labs(title = "Frequency of Life Expectancies across Countries by Continent in 2010", x = "Life Expectancy", y = "Count of Countries", col = "Continent")

Point plots

Point plots display the relationship between a categorical variable and a numeric variable. For instance, the plot below displays a relationship between hospital type and the number of beds at the hospital. Notably, unlike the plots we have been viewing until now, with point plots, we see a point for every observation in the dataset. Because point plots display every observation (rather than aggregating them into other polygons and lines), they are particularly good for seeing outliers in the data. However, with large datasets, this also can mean that points will overlap. Note that the first plot below exhibits overplotting - when the data represented on a plot overlaps, making it difficult to discern one point from the next. There are various tools available to deal with over-plotting. You can reduce the size of points on the plot, increase their transparency, or set them to slightly offset each other (known as adding jitter). We do all three in the second plot below.

#Run this code chunk.

#df %>% ggplot(aes(x = CATEGORICAL_VARIABLE, y = NUMERIC_VARIABLE)) + geom_point()

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = TYPE, y = BEDS)) + 
  geom_point() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  labs(title = "Number of Beds in Hospitals by Type", x = "Type", y = "Number of Beds")

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = TYPE, y = BEDS)) + 
  geom_jitter(size = 0.5, alpha = 0.1) + #Change geom_point to geom_jitter, reduce the size, add transparency for overplotting
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  labs(title = "Number of Beds in Hospitals by Type", x = "Type", y = "Number of Beds")

Notice that think black line in the column for Critical Access hospitals. Based on what we’ve learned of critical access hospitals, can you surmise what number they are clustered at?

(+.25 points on lab if you can respond correctly here.)

Critical access hospitals have many guidelines they have to follow in order to receive funding. In this case, the black line in the column represents the cluster of beds at the value 25. That is because critical access hospitals can't have more than 25 beds in order to adhere to guidelines and receive funding.

Titling a Point Plot

Note how I titled my plot above: “Number of Beds in Hospitals in the US that are Open by Type” Here I am displaying the number of beds by Type. A good formula for titling point plots is as follows:

Number/Measure of [y-axis variable name] in _____ by [x-axis variable name]

The blank should be filled with your unit of observation. The [x-axis variable name] should be your x-label and “Number/Measure of [y-axis variable name]” should be your y-label.

What if I have qualified units of observation?

If this is the case I would encourage you to filter your data as you have been doing elsewhere. For instance if world_health_econ is qualified by country and year, I can filter to the most recent year before creating a point plot of the private share of health spending by continent.

#Run this code chunk.

world_health_econ %>% 
  filter(year == max(year, na.rm = TRUE)) %>%
  ggplot(aes(x = continent, y = priv_share_health_sp)) + 
  geom_point(size = 2, alpha = 0.8) +  
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  labs(title = "Private Share of Health Spending in Countries by Continent in 2010", x = "Continent", y = "Private Share of Health Spending") 

Scatterplots

Scatterplots display the relationship or correlation between two numeric variables. For instance, below we plot the relationship between the BEDS variable and the POPULATION variable in the hospitals dataset.

There are different ways to characterize the correlations between variables in data. We may consider the strength of a correlation, the shape of a correlation, and whether the correlation is positive or negative.

• Two variables in a scatterplot can be said to have a strong correlation when points are clustered closely to a central line or curve. The more scattered throughout the plot the points are, the weaker the correlation.
• Two variables in a scatterplot can be said to have a linear correlation when the scatterplot tends to produce a straight line. This means that the rate of change between two variables is steady. However, when a scatterplot produces a curve, this indicates that the rate of change between two variables is not as constant.
• Two variables in a scatterplot can be said to have a positive correlation when the points move upward from the bottom left corner towards the top right corner of the plot. This means that as values in variable increases, the values in the other variable also increase. When points move downward from the top left corner to the bottom right corner of the plot, we can say that the variables negatively correlation. This means that as values in variable increases, the values in the other variable decreases.

#Run this code chunk.

#ggplot(df, aes(x = NUMERIC_VARIABLE, y = NUMERIC_VARIABLE)) + geom_point()

hospitals %>% 
  filter(STATUS == "OPEN") %>%
  ggplot(aes(x = BEDS, y = POPULATION)) + 
  geom_point(size = 0.5) +
  theme_bw() +
  labs(title = "Relationship between Hospital Beds and Population in the US", x = "Beds", y = "Population") 

Note that this plot has a strong, linear positive correlation. As BEDS increase in this variable, POPULATION also tends to increase.

Titling a Scatterplot

Note how I titled my plot above: “Relationship between Hospital Beds and Population in the US” A good formula for titling scatterplots is as follows:

Relationship between _____ [x-axis variable name] and [y-axis variable name]

The blank should be filled with your unit of observation. The [x-axis variable name] should be your x-label and [y-axis variable name]" should be your y-label.

What if I have qualified units of observation?

If this is the case I would encourage you to filter your plot to one value in one of your qualified variables. For instance if world_health_econ is qualified by country and year, I can filter to the most recent year before creating a point plot of the private share of health spending by continent.

#Run this code chunk.

#ggplot(df, aes(x = NUMERIC_VARIABLE, y = NUMERIC_VARIABLE)) + geom_point()

world_health_econ %>% 
  filter(year == max(year, na.rm = TRUE)) %>%
  ggplot(aes(x = tot_health_sp_pp, y = life_exp, size = pop, col = continent)) + 
  geom_point(shape = 21, stroke = 1) +
  labs(title = "Relationship between Country Total Health Spending Per Person and Life Expectancy", x = "Total Health Spending Per Person", y = "Life Expectancy", size = "Population", col = "Continent") + 
  theme_bw() +
  scale_size_continuous(range = c(1, 10), labels = scales::comma)

Note that this plot has a weaker, curvilinear, positive correlation. As Total Health Spending per Person increases in this variable, Life Expectancy also tends to increase, but the rate of change at which it increases is not constant.

Have you heard the quip “Correlation does not equal causation”? This is particularly important to consider here. In lab 7, we will examine some confounding variables that may be mediating how values appear to correlate in our data. For now, it’s important to note that just because we see a correlation between total health spending and life expectancy does not mean that increasing total health spending in a country causes life expectancy to increase. This is of course a complex issue with lots of other variables involved.

Produce two plots that represent co-variation in your dataset.

You need not include every plot I described above. Be sure to zoom in or out on your data if you have qualified units of observation.

#df %>% ggplot(aes(x = NUMERIC_VARIABLE, col = CATEGORICAL_VARIABLE)) + geom_freqpoly(binwidth = 1)


ipps %>% 
  filter(DRG.Definition %in% (c("313 - Chest Pain", "813 - COAGULATION DISORDERS", "312 - SYNCOPE & COLLAPSE", "652 - KIDNEY TRANSPLANT", "069 - TRANSIENT ISCHEMIA"))) %>%
  ggplot(aes(x = Average.Total.Payments, col = DRG.Definition)) + 
  geom_freqpoly(binwidth = 20000) +
  labs(title = "Frequency of Average Total Payments across Hospitals by Diagnosis", x = "Average Total Payments", y = "Count of Hospitals Diagnoses", col = "Diagnosis Type") + 
  theme_bw() 

What question might this analysis help to address?

Fill your response here. 

Are there any other variables in your dataset that you need to take into consideration before directing this analysis towards answering that question? In other words, do you need to zoom into any specific areas of the dataset (by filtering) in order to appropriately address this question? If so, which? Be sure to adjust your plot above to reflect this.

Fill your response here. 

What else would we need to know to fully address this question? Here you may consider what you know about how this dataset was produced and its limitations.

Fill your response here. 

What insight can you draw from this plot?

Fill your response here. 

#Fill the code for plot 2 here. Add a title and labels to your plots. Be sure to adjust for overplotting.

ipps %>% 
  filter(DRG.Definition %in% (c("313 - Chest Pain", "813 - COAGULATION DISORDERS", "312 - SYNCOPE & COLLAPSE", "652 - KIDNEY TRANSPLANT", "069 - TRANSIENT ISCHEMIA", "038 - EXTRACRANIAL PROCEDURES W CC", "472 - CERVICAL SPINAL FUSION W CC", "551 - MEDICAL BACK PROBLEMS W MCC", "637 - DIABETES W MCC"))) %>%
  ggplot(aes(x = DRG.Definition , y = Average.Total.Payments)) + 
  geom_jitter(size = 0.5, alpha = 0.1) + #Change geom_point to geom_jitter, reduce the size, add transparency for overplotting
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust=1)) +
  labs(title = "Average Total Payments across Hospitals by Diagnosis", x = "Diagnosis Type", y = "Average Total Payments")

NA
NA

What question might this analysis help to address?

Fill your response here. 

Are there any other variables in your dataset that you need to take into consideration before directing this analysis towards answering that question? In other words, do you need to zoom into any specific areas of the dataset (by filtering) in order to appropriately address this question? If so, which? Be sure to adjust your plot above to reflect this.

Fill your response here. 

What else would we need to know to fully address this question? Here you may consider what you know about how this dataset was produced and its limitations.

Fill your response here. 

What insight can you draw from this plot?

Fill your response here. 

Start a shiny app.

At this point in the quarter, we are prepared to start piecing together a dashboard for displaying the data. Shiny Apps have two components - a front end (ui), and a backend (server). On the front end, we will be coding how we want our data displayed. For this week, the only thing you need to do on the front end is fill in your title. On the backend, we will be coding our data analysis.

Fill in your title in the UI. Leave everything else the same.

#Follow instructions to fill code between capitalized sections and then run the code chunk.

ui <- dashboardPage(
  
  #REPLACE 'TITLE HERE' BELOW WITH YOUR OWN TITLE
  dashboardHeader(title = "Medicare Data"),
  #REPLACE 'TITLE HERE' ABOVE WITH YOUR OWN TITLE
  
  dashboardSidebar(
    #inputs will go here eventually but not now. 
  ),
  
  dashboardBody(
      box(plotOutput("plot1")),
      box(plotOutput("plot2"))
  )
)

Replace my plots in the code below with two of the plots that you created in this lab.

#Follow instructions to fill code between capitalized sections and then run the code chunk.

server <- function(input, output, session) {
  
  output$plot1 <- renderPlot({
  #REPLACE BELOW WITH YOUR OWN PLOT
    hospitals %>% 
      filter(STATUS == "OPEN") %>%
      ggplot(aes(x = TYPE)) + 
      geom_bar() +
      labs(title = "Number of Hospitals in the US that are Open by Type", x = "Type", y = "Count of Hospitals") + 
      theme_bw() +
      theme(axis.text.x = element_text(angle = 90, hjust=1)) 
  #REPLACE ABOVE WITH YOUR OWN PLOT
    
  })
  
  output$plot2 <- renderPlot({
  #REPLACE BELOW WITH YOUR OWN PLOT
    hospitals %>% 
      filter(STATUS == "OPEN") %>%
      ggplot(aes(x = TYPE)) + 
      geom_bar() +
      labs(title = "Number of Hospitals in the US that are Open by Type", x = "Type", y = "Count of Hospitals") + 
      theme_bw() +
      theme(axis.text.x = element_text(angle = 90, hjust=1)) 
  #REPLACE ABOVE WITH YOUR OWN PLOT 
  })
  
}

#Run this code chunk to see your Shiny app in action. 

shinyApp(ui, server)

When you knit this notebook, you may get a prompt asking whether you would prefer to render and run this document as Shiny. Click No. We will be leaving this notebook in a HTML document format.